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Abstract. Modified gravity is one of the most promising candidates for explaining 
the current accelerating expansion of the Universe, and even its unification with the 
inflationary epoch. Nevertheless, the wide range of models capable to explain the 
phenomena of dark energy, imposes that current research focuses on a more precise 
study of the possible effects of modified gravity may have on both cosmological and 
local levels. In this paper, we focus on the analysis of a type of modified gravity, the 
so-called f{R, G) gravity and we perform a deep analysis on the stability of important 
cosmological solutions. This not only can help to constrain the form of the gravitational 
action, but also facilitate a better understanding of the behavior of the perturbations in 
this class of higher order theories of gravity, which will lead to a more precise analysis 
of the full spectrum of cosmological perturbations in future. 
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1. Introduction 

In the recent years, modified gravity theories have become one of the most popular candi- 
dates for explaining the current accelerated expansion of the Universe. As it is very well 
known, general relativity (GR) in its standard form can not explain such behavior of 
the Universe expansion without either extra terms in the gravitational Lagrangian (for 
reviews on modified theories of gravity, see Refs. [Ill2j) or exotic fiuid components (see 
Refs. [3]-[5]). Modified gravity theories have been widely studied, and it has been shown 
that they are not only capable to mimic the dark energy epoch, but also the infiationary 
era [6]. Therefore, by the only use of large scale observations (la type supernova, BAO, 
or the cosmic microwave background) which depend uniquely on the evolution history 
of the Universe, the nature and the origin of DE cannot be determined due to the fact 
that identical evolutions for the cosmological background can be explained by a diverse 
number of theories. This is the so called degeneracy problem. It is thus required, in 
order to confirm or discard the validity of these theories, to obtaining solutions that 
can also describe correctly, e.g., the growth factor of scalar perturbations (see Refs. [7J), 
the stability of cosmological solutions against small perturbations and the existence of 
GR-predicted astrophysical objects such as black holes [8]. 

In this sense the simplest, and in fact the most studied, modification of GR is the 
one where the Hilbert-Einstein action is generalized to a general function of the Ricci 
scalar i?, the so-called f{R) gravities [^-[I^. These theories are able to mimic the be- 
havior of the cosmological constant (see for instance Ref. [llj), but also can reproduce 
the entire cosmological history (see Ref. [6]). In addition, they seem to behave quite 
well at local scales, where the GR limit must be recovered [12j, and the existence of GR- 
predicted astrophysical objects such as black holes can be achieved [13j. Nonetheless, 
these candidates have their own shortcomings [lOj and have to pass rigorous theoretical 
and observational scrutiny before they can be accepted as viable theories [12]. An- 
other possible modification of the standard gravitational Lagrangian includes a wider 
number of curvature invariants (i?, Rjj^iyR^^ ^ Rj^xj^ctR^^^^ among others). Within these 
modifications, the so-called Gauss-Bonnet gravity can be included. In these theories 
the gravitational Lagrangian consists of a function /(i?, G) where G holds for the usual 
Gauss-Bonnet invariant. This modified gravity has been also widely studied and it is 
known that can also reproduce any kind of cosmological solution (see Refs. p!4]-[2T]). 
where special attention has been already paid to models able to mimic the ACDM 
model, as well as other important cosmological solutions (see Refs. [T5l[l6]). Finally, 
the cosmological perturbations have been explored within different standard scenarios 
for this class of theories |T7]. 

In this investigation, we are interested in studying the stability behavior of several 
kind of cosmological solutions in the framework of Gauss-Bonnet gravities when sub- 
jected to homogeneous perturbations. Our analysis will therefore exclude anisotropic. 
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i.e. cosmological scalar, perturbations. Homogeneous and isotropic perturbations have 
been historically considered as the mean to determine the stability of different modi- 
fied gravity theories (see for instance Refs. [3l [22l El]). The usual approach in these 
references consisted in perturbing both the Hubble parameter and the matter density 
in order to check the background stability in the time evolution. This way, one can 
determine the stability of the modified Einstein equations when treated as differential 
equations. With respect of the cosmological perturbations, the full anisotropic analysis 
of cosmological perturbations in modified gravity theories is out of the scope of this 
investigation whereas some significant advances have been made in the last years, in 
particular for f{R) theories 0. 

With regard to Gauss-Bonnet gravity theories, a particular subclass has been stud- 
ied in Refs. J^, where the Hilbert-Einstein action plus a function f{G) is considered. 
However, the extension to a more general form for the /(i?, G) gravitational theory is 
a mandatory task that may help to understanding the viability and features of more 
general Lagrangians. These theories may present ghost degrees of freedom in an empty 
anisotropic universe, i.e. the Kasner-type background [20]. Nonetheless, these degrees 
of freedom are absent on Friedmann-Lemaitre- Robertson- Walker (FLRW) backgrounds 
(see also [2^). This is precisely the type of cosmological background that will be con- 
sidered in our investigation. Hence, in the present paper we analyze some important 
cosmological solutions in /(i?, G) gravity, both in a vacuum scenario and with the pres- 
ence of standard perfect fiuids. The stability of cosmological solutions has been studied 
for f{R) gravity in Ref. [22j, as well as for other curvature invariants in Ref. [3]. For 
some particular cases of /(i?, G), the perturbations on the solutions have been analyzed 
in Ref. [23], whereas for Hofava-Lifshitz gravity, the analysis has been performed in 
Ref. [24]. Here we extend the analysis to more general actions, where we shall find 
the stability conditions for different cosmological evolutions in FLRW universes such 
as infiationary epoch and late-time accelerated era as described by the ACDM model. 
Therefore, these analyses can help understanding the viability of cosmological evolution 
provided by this kind of modified gravity, and constraining the viable candidates for the 
underlying gravitational action. 

The paper is organized as follows: in Section H we present the general features of 
the /(i?, G) gravity theories by writing the corresponding modified Einstein equations. 
In Section HI we introduce the evolution equations of perturbations appearing in 
these scenarios once that a FLRW background is assumed. Sections IV and V are 
then devoted to the study of stability around the de Sitter and power-law solutions 
respectively. In the last case, we pay special attention to configurations including 
perfect fiuids such as radiation and dust. Section VI is finally devoted to study the 
stability of the /(i?, G) model able to mimic the ACDM cosmological evolution without 
any cosmological constant. We conclude the paper by giving our conclusions in Section 
VII. An appendix is included at the end of the communication to provide explicitly the 
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2. F{R, G) gravity 

Let us start by writing the most general action for modified Gauss-Bonnet gravity, wliich 
is given by, 

1 



S= (Tx 



J{R,G)+Cr, 



(1) 



where = SttGn^ Gat is the Newton constant and £^ represents the matter Lagrangian. 
The symbol G holds for the Gauss-Bonnet invariant, which is expressed as, 

G = R^- iR.^R^" + R.^xaR'""'" (2) 

Then, by varying expression ([T]) with respect to the metric tensor 5^^^, the modified 
Einstein field equations are obtained ^ 



-2fGR^f"^^R%r - ^fGR''"'''Rpa + 2(V'^VVG)i? ' 2g'^'' (V^ fa) R - A{VpWfG)R''' 



-4(V,VVG)i?^'' + 4(VVG)i?^'^ + ^g''''{VpV,fG)R'^ 
- 4(VpV,/G)i?'^''^'^ - fcR'"' + W/i - g'^'V^fR. (3) 

where V holds for the usual covariant derivative and subindices G and R in f hold for 
derivatives of the gravitational Lagrangian /(i?, G) with respect to those arguments. In 
this investigation, we are interested in studying different cosmological solutions described 
by a fiat FLRW Universe, such that we assume along the paper the metric, 

i=3 

ds" = -dt" + a\t)^{dx'Y (4) 

i=l 

Then, the Hubble parameter H takes its usual definition = a/a, and G and R become 

G = 24{HH^ + //^), R = 6{H + 2H^), (5) 

where the dot holds for the derivative with respect to cosmic time t. By substituting 
expressions ^ and the metric Q into the field equations ([s]), the FLRW equation for 
indices fi = u = yields 

= ^(G/g - / - 24H^fG) + 3iH + H^)fR - SHfn + K^p^. (6) 

Let us assume from now on, that the matter fiuid will be given under the form of 
a perfect fiuid with a constant equation of state = ujprm with the matter energy 
density satisfying the standard continuity equation 

Pm + 3H{l + w)pm = 0. (7) 



On the stability in f{R^ G) gravity 



5 



3. Perturbations of flat FLRW solutions in G) gravity 

Let us now study the homogeneous and isotropic perturbations around a particular 
cosmological solution for the f{R^G) theory described by the action ([T]). Here we 
establish the perturbed equations for a general case, but specific cases will be studied 
in the upcoming sections, specially de Sitter and power law solutions, as well as the 
behavior of ACDM solution in the context of these theories of gravity. 

For this purpose, let us assume a general solution for the cosmological background 
of FLRW metric, which is given by a Hubble parameter H = H^it) that satisfies the 
background equation ^ for a particular f{R^G) model. The evolution of the matter 
energy density can be expressed in terms of this particular solution by solving the 
continuity equation ([T]) yielding 

P™o(t) = Poe-^(^+-™)/^°«^*. (8) 

Since we are interested in studying the perturbations around the solutions H = Ho{t) 
and density given by ([s]), let us consider small deviations from the Hubble parameter 
and the energy density evolution. Hence, 

H{t) = Ho{t) (1 + S{t)) , p^{t) = pmoit) (1 + 5^{t)) . (9) 

where 5{t) and S^it) hold for the isotropic deviation of the background Hubble 
parameter and the matter overdensity respectively. In order to study the behavior 
of these perturbations in the linear regime, we expand the function /(i?, G) in powers 
of Rq and Go evaluated at the solution H = Ho{t)^ 

f{R, G) = f + fl{R - Ro) + fciG - Go) + (10) 

where the superscript refers to the values of /(i?, G) and its derivatives evaluated at 
R = Ro and G = Go- The term includes all the terms proportional to square or 
higher powers of R and G that will be included in the equation, although only the linear 
terms of the induced perturbations are considered. Hence, by introducing expression 
([9]) in the FLRW background equation ^ and using the expansion (10), the equation 
for the perturbation S{t) becomes in the linear approximation, 

C2S{t) + Ci5{t) + CoS{t) = CmSm{t) , (ll) 

where coefficients Co,i,2 and are explicitly written in the Appendix at the end of the 
communication. These coefficients depend explicitly on the /(i?, G) and its derivatives 
evaluated in the background solution. In addition, there is a second perturbed equation 
obtained from the matter continuity equation ([T]) once it is perturbed with expressions 
Thus, 

Sm{t) + 3Ho{t)S{t) = . (12) 
Hence, for a particular FLRW cosmological solution, its stability can be studied in the 



context of /(i?, G) gravity by analyzing and solving the equations (11) and (12). Due to 



the linear character of (11), the solution for 5{t) can be in general split in two branches: 



the first one corresponding to the solution of the homogeneous equation in (11), which 
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reflects the perturbations induced by the chosen particular gravitational Lagrangian. 
The second branch would correspond to the particular solution of that equation, which 
is merely affected by the growth of matter perturbations 5jn{t). Hence, the general 
solution can be written as, 

^ homogeneous {t^ ~l~ ^ particular (j^^ • (-^3) 

In the upcoming sections, we shall consider several cosmological solutions, and their 
stability will be then studied. Nevertheless, let us flrstly do some general considerations. 



By recovering the Hilbert-Einstein /(i?, G) = i?, the stability equation (11) yields, 

- 6H^S{t) = cMt) . (14) 

which can be understood as an algebraic relation between the geometrical and the matter 
perturbations. Hence, in GR the full perturbation around a cosmological solution is 
fully determined by the matter perturbations (or vice versa). In fact, by taking explicit 



expression for in the Appendix and equation (12) it is straightforward to prove that 



m = -ISmit) oc a{tf' (15) 



Nevertheless, the algebraic relation (15) between S{t) and Sm{t) in GR is in general 



absent for higher order theories of gravity. For these theories, the evolution of the 
perturbations is in general determined by a coupled system of ordinary differential 



equations, (11) and (12), where the underlying gravitational Lagrangian plays an 
essential role on the form of the Co,i,2 coefficients as can be seen in the Appendix. 
Let us illustrate the previous comments by considering theories whose Lagrangians can 
be rewritten as follows 

f{R,G) = MG) + MR) . (16) 
Then, the perturbation equation (fTTl) becomes much simpler and some information can 



consequently be extracted. For instance, after assuming that GR should be recovered in 



some limit, which basically means that the higher derivatives of the function (16) should 



become negligible, the coefficient Ci in (11) becomes null, and by diving the equation 
by C2, it yields. 



^^'^ + 3(16 facHoitr + fnn)^^'^ ~ 18(16/°«i/o(t)^ + fnn) ^""^'^ ' ^^^^ 
where cL, = .0/1^ .o S^/^x4 , ^0 n » In order to ensure the stability of a particular solution 
in vacuum in the GR limit and provided that > 0, the denominator in (17) has to 



satisfy, 

16fGGHoit)' + fRR>0. (18) 

In fact, this constraint was also proved in [25] to guarantee the generalized second law 
of Thermodynamics for these theories in de-Sitter scenarios. Note that for /(i?, G) = 
R + /(G), the Lagrangian is restricted to be /qq > to ensure the stability of any 
solution in the GR limit [18j. Nonetheless, this constraint on the second derivative with 
respect to G may produce, according to [19j, instabilities in the matter cosmological 
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perturbations for this kind of models. However, a wider range of functions fi{R) in (16) 
may circumvent the existence of such instabihties. 

4. Stability of De Sitter solutions 

Let us start our study of the stabihty of different cosmological solutions by studying 
some of the simplest cosmological solutions, namely the de-Sitter (dS) solutions, 

Ho{t) = Ho a{t) = aoe^°^ , (19) 

where Hq is constant. According to ([5]), the Ricci and Gauss-Bonnet terms are given in 
this case by Rq = 12Hq and Gq = 24:Hq. Then, inserting the expression for the Hubble 
parameter (19) in equation ^ once vacuum is considered, we get 

^ (Go/G(i?o, Go) - /(Go, Ro)) + SH^ofniRo, Go) = (20) 

Therefore, any f{R^G) function can in principle admit vacuum de Sitter solutions 
provided that the previous algebraic equation has positive roots for Hq. Hence, for some 
particular /(i?, G) models, the current accelerating epoch of the Universe expansion - as 
well as the inflationary epoch - can be explained. Following a different approach, one can 



consider equation (20) as an ordinary differential equation for the /(i?o,Go) function 



so that the corresponding solution would admit any Hq value. Thus, by solving the 



equation (20) in terms of /(i?, G), the following action is obtained, 

/(i?, G) = aG[R- QHl log(G)] (21) 

where a is an arbitrary integration constant. Therefore functions as the one in ( [2T| ) 
admit an inflnite number of vacuum dS solutions Hq. 

Let us now study the stability of de Sitter solutions in vacuum, where the 
perturbation is affected only by the underlying gravitational theory. According to the 



equation for the perturbations (11), it yields, 

- (64i/oVL - \fR + 32i7oV^G + ^Hlfl^ 5{t) = . (22) 

Therefore, the stability of each dS solution depends on the values of the function /(i?, G) 
and its derivatives evaluated in {Rq^Gq}. The general solution of equation ([22]) can be 
easily obtained, and is given by, 

S{t) = Cie^+' + C2e^-' , (23) 

where 



/i± = 9i/oi^° ± V3F0 (_4/o + 75i/2^o). (34) 

where the variable = ^QHq /qq+SHq f^Q+ f^j^ was introduced to lighten the notation. 
The growth of the perturbation will depend both upon the overall sign of the parameters 



fi± in expression (24) and also upon the real or imaginary character of the square root. 



Thus, four different cases can be distinguished: 
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Fq < and < 75HqF^ /A: implies complex solutions and < 0, therefore 

solutions behave as a damped oscillator of decreasing amplitude. Solutions are thus 
stable. 

Fq > and > 75HqF^ /A: implies complex solutions and > 0, therefore 

solutions as an oscillator of increasing amplitude. Solutions are thus non-stable. 

Fq > and ^^^2pQ ^ (0, 1): implies that both ii± are real and /i+ > 0. Solutions 
are thus non-stable. 



Fo < and 



G (0, 1): implies that both /x± are real and /x_ < 0. In this 



case, /i+ > provided that |/^| > |12i7QFo| and //+ < whenever |/^| < |12i7oFo|. 
Therefore the solutions are non-stable and stable respectively. 



According to the Appendix definitions, since F' 



'0 _ 



y||2, then C2 > 0, i.e., the 



positivity of the second order coefficient in the perturbation equation (11), turns out to 



be the necessary (but not sufficient) condition to provide stability of the dS solutions in 
vacuum. 

In order to illustrate the previous calculations, let us consider the function 

/(i?, G) = KrR + K^K'G'^ . (25) 
Here {ki^ K2} are positive coupling constants. For simplicity we take n = 1 and m = 3, 



so that the equation (20) can be solved, and we find the following dS solution, 

H, = ^(^\'\ (26) 



Thus, considering the perturbation solution (24), and introducing the function (|25|) and 



its derivatives evaluated in (|26|), the solution for the perturbation is determined, 

m - 



Ci exp 



X exp 




25^1 ^2 



223/332/,l/3^5/3 



V3^25^ 



Co 



3^+ a/25ki -223/332^ 



o 1/7 2/3 



Kin 



K1K2 



1/6 



-t 



(27) 



Hence, the stability depends on the value of the exponential of the first term in (27), 



which is given by the values of the coupling constants Ki and /^2- In this sense, for 
Ki > 2 (^^)^^^ ^2, the perturbation grows exponentially, and the dS solution becomes 
unstable, otherwise the perturbation turns out to behave as a damped oscillator that 
tends to zero, so that the solution becomes stable. 



5. Stability of power- law solutions 

In this Section we are interested in cosmological solutions of the type, 
a{t) DC r ^ H{t) = - . 



(28) 
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that shall be referred to as power laws. These solutions represent the scale factor 
evolution for standard fluids, such as dust (a = 2/3) or radiation (a = 1/2) dominated 
Universe provided that GR is assumed as the underlying valid gravitational theory. 



For the solutions (28), the Gauss-Bonnet term and the Ricci scalar, given by 



expressions ([7|, take respectively the following form. 



G 



R 



Qa{2a - 1) 



(29) 



t4 ' t2 • 

Let us now explore two different kinds of gravitational Lagrangian /(i?, G) able to 



mimic the power-law solutions described by (28) in the presence of standard fluids. It is 



obvious that one has to solve flrst the background FRLW equation ^ in order to flnd 
the appropriate function /(i?, G). The considered two types of /(i?, G) Lagrangians are 
the following 



5.1. f{R,G) = h{G) + h{R) 

For the sake of simplicity, let us start by considering the subfamily of functions given 
in (16). Then, for this type of Lagrangian the Friedmann equation ([g]) can be split into 
two equations [15j, as 

-2AH'G{h)^a + G{h)^-h 
-3HRif2)j,^ + 3{H + H')if2) ^ 



IR 



:/2 + «VmO 



= 
= 



The first equation in (30) can be written in terms of G, as 



fi = 



which is an Euler equation, whose solution is easily obtained, and yields 
/i(G) = CiG'^+C2G, 



(30) 



(31) 



(32) 



where Ci,2 are integration constants. Note that the second term in (32) is the trivial 



Gauss-Bonnet solution and can be removed, since it becomes a total derivative and does 



not contribute to the field equations. In the same way, the second equation in (30) takes 
the form 

o2 / . N « - 1 ^ . . ^ 2a - 1 



mO 







(33) 



In the presence of a perfect fluid, whose equation of state is given by p^o = '^mPmO, 
from the energy conservation equation ([7| with the class of cosmological solutions given 
in (28), the energy density yields. 



PmO — Pmo{'ttoday)t + _ Pmo{ttoday) 



R 



2 



(34) 



_6«(2tt - 1)_ 

where a 7^ 1/2 has been assumec||} Hence, the equation (j33| for f2{R) gives the solution, 
/2(i?) = hR^+ + k2R^- + , (35) 

J This is a special case to be discussed separately. 
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where ki^2 integration constants and 

3 - tt(7 - 2a) ^1 + tt(10 + tt) 



/5 



A 



4(1 -2«) ^ 4 

-3/t^Pmo(ttodaj/)(2« - 1) 

[6a(2a - 1)]^ [A{A - 1) + - 1) - ^] ' 

3(1 + Wm)0i 



(36) 



Then, the complete function f{R^G) is reconstructed for this class of cosmological so- 



lutions, and can be expressed as the sum of the expressions (32) and (35). 



Concerning the stability of this kind of solutions, we have to evaluate the function 



/(i?, G) and its derivatives in (28), and solve the perturbation equations (11) and (12). 



However, analytical solutions for equation (11) are in general difficult to be found. Only 



numerical solutions can be obtained by fixing the values of the free parameters of the 
theory, i.e., the coupling constants. In order to circumvent this intrinsic difficulty, let us 
consider the following particular cases: the radiation and dust matter evolutions, where 
a = 2/3 and a = 1/2 respectively, and a > 1, which gives an accelerating expansion. 



5.1.1. Dust dominated Universe^ a = 2/3: In the case of an expansion of the type 
of dust matter, given by a = 2/3, the Gauss-Bonnet term G < according to (29), 
so that the gravitational Lagrangian (32) becomes complex, which may be interpreted 
as a non-physical case. In order to avoid this scenario, we set Ci = 0, and thus the 
Lagrangian turns out /(i?, G) = /2(i?), given by (35) . For solving the differential 
equation (11), we use numerical methods. In order to illustrate the richness of this 
case, we considered specific values for the coupling constants ki = k2 = O.lp^o with 
the appropriate dimensions. In figure [l| the evolution of 6{z) and Sm{z) are plotted as 
functions of redshift z, and assuming different initial conditions for both perturbations. 
We also explore the effects of the variation of the initial condition for the first derivative 
S\z) in figure [2| We can see that regardless the initial conditions, the perturbation S 
oscillates tending to zero at z = and similarly for 5^ in figures [T] and [2] Moreover, 
the case of a radiation-like fiuid = 1/3 is explored in figures [3] and [ij where the 
perturbations increase at small values of the redshift z = in comparison with the 
pressureless case. While in figure [5] the values S{z = 0) and Sm{z = 0) are plotted versus 
the initial conditions assumed at z = 1000. 



5.1.2. Late-time acceleration^ a > 1: Another interesting example is given by a power- 
law cosmological solution with a > 1. In GR, this kind of solution is provided by the 
presence of a perfect fiuid whose EoS parameter is given hy w < —1/3 and consequently 
compatible with the currently observed accelerating expansion. Coming back to /(i?, G) 
gravity theories, we are interested in studying its stability when dust matter = 
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Figure 1: The stability of the solutions for the model f{R,G) = fi{G) + f2{R) is shown. Here we 
assume a = 2/3 and Wm = 0, and the action reduces to f{R, G) = f2{R)' The perturbations S{z) and 
6m{z) are represented versus the redshift z. At the top figures, the initial condition dmo at z = 1000 
is varied while So = 0.001 and Sq = are fixed. At the bottom, (5^(1000) = 0.001 is fixed while 
the evolution of 6{z) and Sm{z) is shown as a function of Sq. As shown, independently of the initial 
conditions, the perturbations behave as a damped oscillator, turning out very small at z = 0. 




Figure 2: As Fig. [l| the stability for the solutions described by a = 2/3 and Wm = is analyzed. 
Here the perturbations S{z) and Sm{z) versus the redshift z are represented and the initial conditions 
^0 = ^mo = 0.001 are fixed, while 5q is varied along the range { — 10~^, 10~^}. Unlike Fig. [l] the values 
of the initial condition Sq may produce a large amplitude in the oscillations of the perturbations along 
the evolution, leading to large effects at z = 0. 



is included. It can be shown that the equation (11) exhibits a pole, whose position 
in relation with the redshift depends upon the coupling constants values. For small 
redshifts the perturbation oscillates close to z = 0, where the dark energy epoch is 
expected. Hence, it seems that in the presence of dust matter (baryons and cold 
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Figure 3: The stability of the power law solution with a = 2/3 and Wm = 1/3 is studied for the 
model (35). As in Fig. [l] at the top figures the initial condition 5mO at z = 1000 is varied whereas 



^0 = 10~^ and = are fixed. At the bottom, ^^(1000) = 10~^ is fixed while the evolution of S{z) 
and 6m{z) is shown as a function of ^o- Here, the perturbations for non null initial conditions grow 
along the redshift evolution, producing instabilities at z = 0. 



O.ftftOO 



, ft.OOOS 



0.0010 



-0.0010 




loga + 1) 3 



log(z + 1) 



Figure 4: As in Fig. [sj the stability for the solutions given by a = 2/3 and Wm = 1/3 is analyzed. 
However, here the perturbations 5{z) and Sm{z) are represented versus the redshift z with the initial 
conditions = ^mo = 0.001, whereas Sq is varied along the range { — 10~^, 10~^}. In this case, the 
perturbations grow very fast for non null initial conditions on (5q, producing large instabilities at z = 0. 



dark matter), this kind of solutions becomes unstable for large redshifts, when the 
radiation/matter dominated epochs occur, and naturally the power-law with a > 1 is 
not valid anymore. However the solution tuns out to be stable when we approach to the 
current epoch, and late-time acceleration is well reproduced. For illustrative purposes, 
we have considered fci = /c2 = 0.1 x po and Ci ^ PmoiUoday) whose perturbations are 
depicted in Fig.lGl In this case the contribution coming from the Gauss-Bonnet part (32) 
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Figure 5: This represents a summary of the stabihty of f{R, G) = fi{G) + f2{R) for the power law 
solutions in the dust case a = 2/3: Evolution of cosmological and dust perturbations S and Sm for 
= (left panel) and Wm = 1/3 (right panel) respectively. Initial conditions {So^Smo} are imposed 
at redshift z = 1000 and varied in the range {-10"^, 10"^}, while S\z = 1000) = is fixed. In both 
figures, the hypersurfaces with larger slope correspond to Sm{z = 0), while the other hypersurfaces 
corresponds to 5{z = 0), which are not constant but with a much more little growth than Sm- In both 
cases, the evolution of the final values for the perturbations follow a linear relation with the initial 
conditions, as natural since the perturbations equations are studied in the linear regime. 




Figure 6: Stability of f{R, G) = fi{G)^f2{R) for a = 2: The evolution of the perturbations S{z) and 
Sm{z) versus the redshift z, and the initial values for 5mO at the top, and Sq at the bottom. Here we 
consider a = 2 and initial conditions imposed close to the pole z ~ Zpoie where Sq = 0.001 and Sq = 
are fixed for the top figures, and S^o = 0.001 for the bottom figures. In both cases, we have assumed 
Wm = 0. Different values of the coupling constants change the position of the pole Zpoie- In all the 
cases the perturbation S{z) behaves smoothly close to 2; = while it obviously becomes very large for 
z close to the pole. In the same way 5m{z) tends to small values close to z = 0. 



On the stability in f{R^ G) gravity 



14 




Figure 7: Evolution of cosmological and dust perturbations 5 and 5m for a = 2 with the presence 
of a pressureless fimdwm = 0. Initial conditions {(^c^mo} are imposed at redshift z = and varied 
in the range { — 10~^, 10~^}, while (^'(0) = is fixed. As above, the hypersurface with larger slope 
corresponds to 5m{z = 2.08), while the other hypersurface corresponds to S{z = 2.08). The final values 
are calculated at z = 2.08, very close to the pole shown in fig. [6) The perturbations do not become 
larger than the initial values assumed before the divergence. 



has to be included since G > and the action is therefore real. As in the cases above, 
different initial conditions for both S and were considered, but now at z = since we 
are interested in studying the past evolution of this solution with a > 1, as it reproduces 
late-time acceleration at the current epoch. As shown in Fig.|6} the perturbations diverge 
for large redshifts, so the solution is stable around z = but introduces large instabilities 
for large redshifts, where naturally the presence of a pressureless fluid (dust) should 
dominate. Moreover, the values of 6{z = 0) and 5rn(^ = 0) are evaluated depending on 
the initial conditions S{z = 2.08) and Sm{z = 2.08) in figure [?} 



5.1.3. Radiation dominated Universe^ a = 1/2; In the case for a = 1/2, the equation 
(33) has to be solved separately from the general case, which for a = 1/2, becomes, 

R'if2)RR-\Rif2)R = , (37) 

whose general solution is given by, 

f{R) = CrR'^'^ + C2 . (38) 

Nevertheless, note that this is not the only solution for the case a = 1/2 since this 
value for the power-law exponent implies straightforwardly i? = 0. Thus, any function 
f2{R) that accomplish limi^^o/2(^) ^ R^ with n > 0, will satisfy the equation (37). 
However, by analyzing the Friedmann equation for f2{R) in (30), in terms of the cosmic 
time instead of the Ricci scalar, one can show that no gravitational action other than 
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Hilbert-Einstein (n = 1) will satisfy this equation. Let us prove the previous statement 
by considering an hypothetical action given by f2{R) oc i?^, the equation (30) yields, 

n-l 



l-2a + n{2n + a - 3) 
t^ 



a{2a- 1) 

t2 



(39) 



Then, if a = 1/2 is assumed, the first term in the l.h.s. of (39) vanishes, whilst the 
second is never zero unless p^o = (vacuum). Hence, the only possible solution is given 
for n = 1, where the f2{R) action becomes the Hilbert-Einstein action, which imposes 
that = 1/3 as natural for a radiation type expansion. 



Hence for a = 1/2, the gravitational action needs to be either /(i?, G) = i? + /i(G), 
or /(i?, G) = f{G). In the first case, for fi{G) given by expression (32), but since G < 



according to ([29|), the action ([32|) would become complex, and non physical, such that 
/(i?, G) = R. For the second possibility, f{R^G) = /(G), the action only depends on 



the Gauss-Bonnet term, and the FLRW equations (30) can be rewritten as 

1 



G fcG + 



a 



-Gf, 



a 



G 



1 







4 4 4 

while the expression for the energy density of a perfect fluid Pm = Wmpm 

3(l+tum)a 



, , , , 3(1 + Wm) 

Pm(t) = PmO(ttoday)t ^ 



PmO 



G 



24a3(a - 1) 



Hence, the equation (|40j) can be easily solved, 
/(G) = CiG'-^ + G2G + /3G^ , 



(40) 
is given by, 

(41) 
(42) 



where 



p = 



A 



f^'^ PmO {t today) 



a 



4 [2ia^{a - 1)]^ [{A - 1)A + ^{A - 1)] ' 

3(1 + Wm)a 



(43) 



However, as above the Gauss-Bonnet term G < since a < 1 and the gravitational 
action becomes complex, which in principle lacks any physical sense. Hence, in the 
case oi a = 1/2, the gravitational action reduces to the Hilbert-Einstein action with 
the presence of a radiation-like fiuid. This analysis provides an interesting consequence, 
since the universe evolution crosses a radiation dominated epoch, if the gravitational 
action is of the type of /(i?, G) = /i(G) + /2(^), the extra terms in the action should 
be negligible during such epoch, and the action must approach Hilbert-Einstein action. 

5.2. f{R,G) = fiR^G^ 

Let us assume now a gravitational Lagrangian of the form 

f{R,G) = fiR^G^ (44) 
where /i is a dimensionful constant and /3 and 7 are constant exponents an let us study 



the power-law solutions as the ones given at (28) for two standard fiuids, radiation an 
dust matter. 
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5.2.1. Radiation dominated era, a = 1/2: For f{R,G) models as given in expression 



(44), it is straightforward to prove that in the absence of any fluid and provided that 



P > 1 and any even value of 7 (since G < 0, this constraint gets rid of imaginary 
terms), they can mimic radiation- like scale factor evolution, i.e. a = 1/2. This can 
be performed by solving equation ^ in the absence of fluid sources and considering a 
function like (44). For this kind of models, equation (11) becomes 



t^S"{t) - It (4/5 + 87-9) S'{t) - ^ (4/5 



2 V ' ■ ' / V / 2 
This is an Euler-like equation that yields to the solutions 

5{t) = C+ r+ + C7_t"- 
where C± are arbitrary constants and 

«± = /3 + 1 (-7 + 87 ± I87 + 4/5 - 3|) 



87 - 5) S{t) = 



(45) 



(46) 



(47) 



In order to understand the stability of the perturbation given by expression ([46|), one 
should study the sign of a± exponents: provided that the requirement /3 > 1 is satisfied, 
it can be proved that is always negative. Concerning this exponent is negative 
provided that 7 < 1/8(5 — 4/3) and 7 even number as explained above. Otherwise 
would be positive and the perturbation unstable. 

If one now considers a radiation fiuid characterized by both uj = 1/3 and a = 1/2 
it is possible to show that Lagrangians such as (44) cannot satisfy equation ^ for any 
value of their parameters. If we relax now the constraint on a permitting a ^ 1/2, then 
the combination /3 + 2^ has to be negative in order to guarantee a > . It can be 
shown that in this case, equation ([g]) is not satisfied for any parameter combinatioi][|} 
Therefore, none of these two last cases can be accomplished by Lagrangians such as 



(44) 



5.2.2. Dust dominated era, a = 2/3: For this case, we have considered three different 



scenarios: In the absence of matter, models given by (44) can hold a power-law scale 
factor with a = 2/3. In this case, /i can take in principle any value whereas 7 (which 
has to be even as in the previous case) and /5 must be related as follows 



7 



= 7± = ^(l3 + 6/5±^12r 



180/5 + 324/52 



(48) 



in order to satisfy the background equation With this requirement, the perturbation 



equation (11) becomes a Euler-like equation 

1 

6 

1 



t^S{t) 



+ 



6 



17 + 18/5 ± ^121 + 36/5(-5 + 9/5) t5{t) 
5 - 18/5 T ^121 + 36/5(-5 + 9/5)1 S{t) = 



(49) 



In fact, the required /i would be imaginary. 
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Figure 8: Values today {z = 0) for cosmological and matter (dust) perturbations, 5 and 5m 
respectively for /3 = 5 (left panel) and /3 = — 1 (right panel) in expression (53) and i^o(^) = 2/3t 
in (12). Initial conditions were imposed at redshift z = 1000 ranging in the interval (—0.001, +0.001) 
for both 8{z = 1000) and 5m{z = 1000). In both panels, the hyperplanes with larger slope correspond 
to 6m evaluated today. The other hyperplanes which seem of constant value in the 3-dimensional 
representation correspond to S today. In both panels can be seen how Sm achieves today values bigger 
- in absolute value - than the initial conditions. On the other hand, S remains with amplitudes today 
smaller than the initial values. 



whose solutions are 



t 



(50) 



where the second exponent is always a real number regardless the /3 value. Moreover, 



it can be proved that for i.e. 7+ choice, the second power-law in (50) possesses 



a positive exponent, whereas holds both power-law solutions with negative expo- 

nents. Therefore, the solutions provided by the 7_ choice would be stable whereas the 
ones by 7+ would not be. 



The second case consists in considering the presence of dust matter {uu = 0) and a 
power-law scale factor with exponent a = 2/3. In that case 7 and /i depend on /3 as 
follows 

^ = ^ ' ^ = • ^^^^ 

From the last expression, it is straightforward to conclude that one of the two following 

constraints must be imposed to guarantee fi to be positive and real: 

5 5 
/5 = 4n + 1 and /3 > - ; /3 = 4n + 3 and /3 < - with neZ (52) 

9 9 

For this case, the equation (11) is again Euler-type with a source term proportional to 
(5m (t). Thus 



t^S{t) + 3tS{t) 



1 



12 



45(-l + /3) 5 + 45/5 



m = 



2(-5 + 9/3) 
9[1 + (8 - 9/5)/5] 

(53) 



SM{t) 
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Figure 9: Stability of f{R^G) given by (44) in the dust {uj = 0) case: Evolution of cosmological and 
dust perturbations, S and 5m respectively for /3 = 5 in expression (53). Different initial conditions at 
redshift z = 1000 were imposed. The two upper panels represent the evolution in redshift for S (left) and 
Sm (right) for fixed Sm{z = 1000) = 10"^ whereas 5{z = 1000) varied from -0.001 to 0.001. One can 
see how the evolutions of both perturbations acquire decreasing amplitudes when approaching today. 
At the lower figures, it is depicted again S (left) and 5m (right) this time for fixed S{z = 1000) = 10~^ 
whereas Sm{z = 1000) varied from —0.001 to 0.001. On the left panel it can be seen how S tends to 
null amplitude today regardless the values for 5m{z = 1000). On the contrary, the values of Sm (right 
panel) tend to decrease in amplitude by depending on its initial value. 



In order to illustrate the rich phenomenology of this case, we have considered /3 = 5 
for the first constraint[||] and /5 = — 1 for the second according to (52). For these two 
cases, in Figure [8] we have plotted the values for 5 and today for a wide range of initial 
conditions. On the one hand, for /5 = 5, the homogeneous Euler-like associated equation 
to (53) would present complex conjugate exponents. The presence of the matter term 
induces solutions as the ones presented in Figures ^ In this figure, the perturbations 
amplitudes for 6 seem to decrease whereas those for 5^ increase but remaining all of 
them in the linear regime. On the other hand, for /3 = — 1 the homogeneous associated 
equation to (53) would present two distinct negative real exponents. The inclusion of 
the matter term leads to the perturbations amplitudes to decrease in absolute value as 
presented in Figure 10 both for matter and Hubble parameter perturbations. 

Finally, we decided to study a dust fluid characterized by both uj = but a 2/3. 
It is possible to show that Lagrangians such as (44) can satisfy equation ^ provided 

II Note that /3 = 1 will imply 7 = 0, i.e., the usual EH Lagrangian would be recovered. 
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Figure 10: Stability of f{R, G) given by (44 ) in the dust (a; = 0) case: Evolution of cosmological and 
dust perturbations, 8 and 8^ respectively for /3 = —1 in expression ([53|. Different initial conditions at 
redshift z = 1000 were imposed. The two upper panels represent the evolution in redshift for S (left) 
and Sm (right) for fixed Sm{z = 1000) = 10"^ whereas S{z = 1000) varied from -0.001 to 0.001. One 
can see how the S perturbations tend to decrease its amplitudes when approaching today. Concerning 
6m, with initial amplitude of 10~^, acquired final amplitudes ranging from to 10~^. At the lower 
figures, it is depicted again S (left) and Sm (right) this time for fixed S{z = 1000) = 10~^ whereas 
Sm{z = 1000) varied from —0.001 to 0.001. On the left panel it can be seen that S tends today to 
amplitudes ranging between zero and —5 • 10~^, i.e., the amplitudes decrease for all the Sm initial 
conditions. Concerning the evolution for Sm (right panel), the amplitudes today are smaller than the 
corresponding initial amplitudes. 



that 



a 



3 (/3 + 27) 



(54) 



and /X a certain combination of {/3^^^ PmiUoday)} that guarantees this parameter to be 
both real and positive. Therefore, models like (44) may host dust fluid (u = 0) with a 
power-law evolution exponent a 7^ 2/3 for a suitable choice of parameters. In order to 
illustrate this case, we consider one example with /3 = — 1 and 7 = 3/2 and consequently 
/i = 75 / 256 ^3/ 2/^^p^ (ttoday) ' According to expression (54), the obtained value for the 
a exponent is 4/3. Figure 11 represents 6 and 5^ evaluated today for a wide range of 
initial conditions for these two quantities. On the other hand. Figure 12 represents the 
redshift evolution also for S and 5^ by fixing different initial conditions. There, it can 
be seen how whereas the Hubble parameter perturbation remains small and decreases in 
amplitude, the matter perturbation grows in amplitude whereas remaining in the linear 
regime. 
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Figure 11: Values today for cosmological and matter (dust) perturbations, 5 and 5m respectively, 
for the power-law scale factor with exponent a = 4/3: Parameters of the gravitational Lagrangian 
were chosen to be /3 = —1, 7 = 3/2 and /i = 75 / 256 y^3 /2i<i^ {Uoday)- The hyperplane with larger 
slope corresponds to 5m evaluated today whose values range from —0.03 to 0.03, i.e., 30 times bigger 
than initial amplitudes and therefore showing the instability of the 5m evolution. With regard to 5^ 
the hyperplane for 5 today acquires amplitudes ranging from —3 • 10~^ to 3 • 10~^, i.e., three times the 
initial value and consequently showing as well the instability of the 5 evolution. 



6. Stability of f{R,G) mimicking ACDM solution 

We analize in this section the stability of the f{R^G) function mimicking ACDM 
cosmological evolution without any cosmological constant term. This model was 
originally presented in Following [15] and its key features will be described below. The 
cosmological speed-up effect of the cosmological constant in the Concordance model is 
precisely replaced by the modification introduced by /(i?, G) with respect to the usual 
EH Lagrangian. 

The scale factor solution within GR with dust and cosmological constant is given 

by: 

''(*'-(|)'''-*"'(^^°') (^^' 

where Hq holds for Hubble parameter today and Qrn,A hold respectively for dust and 
cosmological constant fractional densities todajj^ According to [15j, the Gauss-Bonnet 
contribution to the gravitational action that is able to mimic ACDM evolution is given 
by 

f{G) = 6 + ^ C(G) + Hie , (56) 

% For illustrative purposes we shall consider Qm = 0.27 and = 0.73. 
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Figure 12: Stability of f{R,G) given by (|44| in the case of dust {u = 0) and scale factor with 
power-law exponent a = 4/3: Evolution in redshift of cosmological and dust perturbations, S and Sm 
respectively. As in Figure [TT] parameters of the gravitational Lagrangian were chosen to be /3 = — 1, 



3/2 and /i = 75/256^s/3/2hi^ pm{ttoday)- Different initial conditions at redshift 



7 



1000 were 

imposed. The two upper panels represent the evolution in redshift for S (left) and 6^ (right) for 
fixed Sm{z = 1000) = 10"^ whereas S{z = 1000) varied from -0.001 to 0.001. One can see how the 
S perturbations tend to reduce its amplitudes when approaching today showing the stability of this 
case. Concerning with initial amplitude of 10~^, acquired final amplitudes ranging from 10~^ to 
4 • 10~^,i.e., increasing its amplitude by a factor 4. At the lower figures, it is depicted again S (left) and 
Sm (right) this time for fixed 5{z = 1000) = 10"^ whereas Sm{z = 1000) varied from -0.001 to 0.001. 
On the left panel it can be seen how 5 today acquires values ranging from to 5 • 10~^, i.e., decreasing 
its amplitude (initially 10~^) for the 5m initial condition. Concerning the evolution for Sm (right panel), 
the amplitudes today range from (for Sm{z = 1000) = 10"^) to -4-10-^ (for Sm{z = 1000) = -10"^). 



where 



CiG) 



6/ 



(57) 



with I = Hi^ pQa{t = ttoday) ^/3 and G is given by expression ([5]), while according to [15j, 
{6^, i9} are 



e 



^0 L 



9 



[e + D 



d = I 



-K 



(58) 



where e is a constant]^ as well. Therefore, the full gravitational Lagrangian is expressed 
as follows [15l 



/(i?, G) = R+-{e aCf + d aO) + Hie) 



(59) 



+ In the original reference, authors used 5 symbol for this constant. In order to avoid confusion, we 
have preferred to use e symbol. 
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Now that we have revised the form of the /(i?, G) Lagrangian able to mimic the 



background solution (55), we study the system made of by the equations (11) and (12). 



The stability of this model for several different initial conditions can be established has 



been represented in Figs. 13, 14 and 15, Both signs in expression (57) were considered 



leading to two different analyses: For the positive branch in expression (57), the obtained 



solutions were oscillatory with decreasing amplitude for all the studied initial conditions. 
According to the numerical results depicted in Figures flS] (left panel) and 14, Smif) 



attains bigger amplitudes today (in absolute value) than 5{t) in this branch. On the 



other hand, for the negative branch in expression (57), the 5 decrease in amplitude 
whereas 5^ solutions increase for all the studied initial conditions. According to Figure 



13) (right panel), 5^ increases in absolute value faster than 5. Fig. 15 illustrates how 5 
amplitude tends to decay whereas 5^ increases. 





Figure 13: Values today for cosmological and dust perturbations, 5 and 5m respectively, for the 



f{R^G) model (59) mimicking ACDM. Positive (negative) sign in expression (57) were represented in 
left (right) panels. Initial conditions were imposed at redshift z = 1000 ranking from (—0.001, +0.001) 
for both S{z = 1000) and Sm{z = 1000). For the positive branch (left panel) the maximum amplitude 
of 5 today is ±2 • 10~^. Thus, the final amplitude is 20 times smaller than the initial ones guaranteeing 
that S remains in the linear regime and the perturbations remain small. The final value for 6 turns out 
to be independent of the initial conditions for S and does depend solely upon the initial conditions of 
Sjn- Concerning (5^, its final amplitude for Sm acquires maximum-minimum values of ±2 • 10~^, i.e., 
twice times the initial amplitude. These values depend both upon initial conditions for S and S^- For 
the negative branch (right panel) the maximum amplitude of S today is ±2 • 10~^. The value for 
today depends both upon the initial conditions fixed for S and 5^- The maximum amplitude for 5m is 
±0.04, i.e, 40 times bigger than the initial amplitude. 
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Figure 14: Stability of f{R^G) given by (59) for positive sign in expression (57): Evolution in 
redshift of cosmological and dust perturbations, S and 6m respectively. Different initial conditions at 
redshift z = 1000 were imposed. The two upper panels represent the evolution in redshift for 6 (left) 
and Sm (right) for fixed 6m{z = 1000) = 10"^ whereas 6{z = 1000) varied from -0.001 to 0.001. On 
the left panel, one can see that regardless the initial conditions for S this quantity today approaches 
null amplitude with decreasing oscillating amplitude. On the right panel, Sm follows the same behavior 
decaying from 10~^ to null amplitude today. At the lower figures, 6 (left) and Sm (right) are depicted 
again this time for fixed S{z = 1000) = 10"^ whereas Sm{z = 1000) varied from -0.001 to 0.001. The 
oscillatory character of the upper figures appears again. In fact, in both panels of the lower figures, 6 
tends to zero regardless the initial value of Sm as well as does Sm today. 
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Figure 15: Stability of f{R^G) given by (59) for negative sign in expression (57): Evolution in 
redshift of cosmological and dust perturbations, S and Sm respectively. Different initial conditions at 
redshift z = 1000 were imposed. The two upper panels represent the evolution in redshift for S (left) 
and Sm (right) for fixed Sm{z = 1000) = 10"^ whereas S{z = 1000) varied from -0.001 to 0.001. On 
the left panel, one can see that regardless the initial conditions for S this quantity today approaches 
null amplitude. On the right panel, 5m instead grows to bigger amplitudes achieving today values 
between 2 — 4 • 10~^. At the lower figures, S (left) and Sm (right) are depicted again this time for fixed 
5{z = 1000) = 10"^ whereas 5m{z = 1000) varied from -0.001 to 0.001. 5 tends to zero amplitude 
regardless the initial value of Sm- Nonetheless, Sm increases its final amplitude today ranging from 
—4 • 10~^ to 2 • 10~^ with regard to the initial amplitude. 



7. Conclusions 

In this work we have extended the study of /(i?, G) modified gravity theories to provide 
an accurate description of the instabilities that these theories may present. For such 
scenarios, the linearized perturbed equations have been derived. They were obtained 
once the modified Einstein equations are implemented with perturbations in both the 
Hubble parameter and matter density. The resulting coefficients have been explicitly 
presented for the first time in the existing literature and a strong dependence on the 
chosen /(i?, G) model was observed. 

We have studied three of the most important cosmological solutions in the standard 
cosmological concordance model around a spatially fiat FLRW background: de Sitter 
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expansion, power laws and the scale factor solution as provided for ACDM model. The 
required /(i?, G) gravitational Lagrangians to provide de Sitter and power-law solutions 
are explicitly determined whereas the model mimicking ACDM is revised . Concerning 
the perturbations, in the first case, we have found that a certain combination in the 
sign of some derivatives of the gravitational Lagrangian evaluated in the cosmological 
background is required to ensure the stability of the solution. This may provide a way 
to understand the end of an inflationary era produced by one of these de Sitter points 
(which appear as natural solutions of /(i?, G) gravity), as well as for understanding the 
evolution of dark energy epoch. 

With respect to the power-law solutions, the complexity of the coeflicients for the 
perturbed equations led us to study two representative models: either a sum of functions 
of scalar curvature and Gauss-Bonnet term or powered products of these two scalars. 
These models may encapsulate the main features for cosmologically viable /(i?, G) 
theories in some asymptotic regimes. Their rich phenomenology was summarized in 
Table [H 

Besides, the usual power-law behavior for the scale factor in general relativity when 
either dust or radiation are the dominant components in the fluid sector, was found to 
be mimicked, even in the absence of such fluids, by appropriate choices in the pa- 
rameter space of these models. The radiation dominated evolution, i.e. a{t) oc t^/^, 
deserves special attention in the case of actions of the type f{R^G) = f{G) + /(i?). 
For this case, we have shown that the background evolution requires the action reduces 
to Hilbert-Einstein action, a result that was also naturally extended to f{R) theories. 
Consequently, the radiation dominated era can not be described in principle for this 
kind of Lagrangians, or at least, effects of such terms should be negligible during that 
epoch. Concerning the evolution of the perturbations, we have shown that for models 
of the form /(i?, G) = fi{G) + /2(^), the perturbations oscillate tending to zero at 
z = provided that the initial condition 6q is assumed to be null, and the background 
evolution (a) coincides with the EoS parameter of the perfect fluid. Also note that 
independently of the model, the flnal values has a direct linear relation with the initial 
ones imposed. Note also that in the case a = 2 (accelerating expansion), the perturbed 
equation presents a pole at a particular redshift, while the perturbations remain small 
up to the divergence. 

Finally, the case of the /(i?, G) model able to mimic the ACDM scale factor led us 
as well to relevant conclusions. According to previous literature, this model presents two 
branches. Our analysis showed the decreasing and oscillatory character of perturbations 
in the positive branch regardless the initial conditions. Therefore, this solution may 
be considered as stable with respect to small perturbations. On the other hand, the 
negative branch showed that matter perturbations are not constrained in amplitude 
even if the perturbations for the Hubble parameter approach zero for a wide choice of 
initial conditions. 
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Hence, the wide study on both cosmological solutions and stabihty of homogeneous 
perturbations for in /(i?, G) gravity carried out in this paper may help to a better 
understanding for such higher order theories of gravity. We have shown that 
gravitational action plays a very important role in the stability of the solutions depending 
both upon the form of the /(i?, G) theory and the parameters of the model, which means 
a great difference with respect to general relativity. The search for stability in widely 
accepted cosmological solutions helps to constrain the parameters space of the /(i?, G) 
models that may be viable and therefore may deserve further study by future analyses 
of the full spectrum of cosmological perturbations in modified gravity. 



Model /(i?,G) = /i(G) + /2(i?) 


a{t) oc r 


EoS 


Stability 


2 
3 


Wm = Q 
Wm = 1/3 


Stable with = 0. 5^ grows very fast when 7^ 0. 
Unstable: 5 and 5m grow with z. 


2 


Wm = 


Existence of a pole Zp^e fixed at large redshift. Stable close to z = 0. 



Model f{R, G) = ^iR^G-^ 


Configuration 


a{t) oc t« 


EoS 


Stability 


Vacuum 


a = i 

a = l 




/5 > 1, 7 one or even number and 7 < | (5 — 4/3) 
Any /3 7 one or even number and 7 = 7_ given by eqn. (48) 


Fluid 


a = 1 
a = l 

a = l 


- = l 

- = l 

uj = 


Impossible to satisfy. 
Eqns. (47) and (48) must be satisfied, 
5m decreases and 5 increases. 
Impossible to satisfy. 
= — 1,7=|,5 tends to zero and 5^ increases. 



Table 1: Summarization of the stability of power-law solutions for f{R^ G) models given by expressions 
(16) and (44). In the cases where no fluids are considered (Vacuum) it can be seen that stability of 
the cosmological solutions can be achieved for appropriate choice in the parameters space. Once either 
dust of radiation fluids are considered, one sees that usual power-law evolution is on the one hand not 
feasible for radiation and, on the other hand, the stability depends on the initial conditions for dust. By 
the a particular action for this kind of cosmological solutions, modified gravity may contribute during 
the matter/radiation dominated eras, but also reproduce dark energy epoch(a = 2). 
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Appendix 



In this Appendix we present explicitly the coefficients of equation (11). 



C2 



- ISHq |l536i7o^o/3G - IrrHo - ^SH^ /qq + 8Ho{3f^QQ + 2Hof^Q) 



Irr + ^H^Urg + + ^IrrgHq) - SSAH^Hof^Q 



:Co = - 



18432H','Hof^a + 3i^oVL + 192^0 [Igg " ^^M^Irgg + ^IsgHo 

4i/o [t/L + 6/3^^ + sM^fac + Wrrg + WrggHo)] } 
>^ + 3i7o(7/L + 8Ho{^f^G + ^fsR + Q/rrgHo))] - m8H',HofsG 
34567^0 ^o(/^GG + IsgHq) - 288i/oi/o(/GG' + ^Irrg + ^IrggHq) 



+ 12i/o' \ Irr 



Hi 



Cm = I^^Pmoit) . 



(60) 



